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Abstract. The crystal structure of solid oxygen at low temperatures and at pressures 
up to 7 GPa is studied by theoretical calculations. In the calculations, the adiabatic 
potential of the crystal is approximated by a superposition of pair-potentials between 
oxygen molecules calculated by an ab-initio method. The monoclinic a structure is 
stable up to 6 GPa and calculated lattice parameters agree well with experiments. The 
origin of a distortion and that of an anisotropic lattice compressibility of the basal plane 
of (X-O2 are clearly demonstrated. In the pressure range from 6 to 7 GPa, two kinds 
of structures are proposed by X-ray diffraction experiments: the a and orthorhombic 
5 structures. It is found that the energy difference between these structures becomes 
very small in this pressure range. The relation between this trend and the incompatible 
results of X-ray diffraction experiments is discussed. 



PACS numbers: 71.15.Nc, 61.50.Ks, 75.50.-y, 61.66.-f 
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1. Introduction 

At low temperatures or under pressures, molecular oxygen is solidified by weak 
intermolecular interactions. At zero pressure, oxygen transforms to monoclinic a-0 2 
through 7- and P-O2 by cooling. In the a phase, oxygen molecules condense with its 
molecular axis perpendicular to the basal plane of the monoclinic lattice. As shown in 
figure 1(a), experiments demonstrated that the unit cell includes two oxygen molecules 
and the structure belongs to C2/m space group pQ. The ground electronic state of the 
oxygen molecule is spin triplet state. The molecular spins, which are perpendicular to 
the molecular axis, order antiferromagnetically on the basal plane of OJ-O2 with easy 
axis parallel to the 6-axis. The basal plane is illustrated in figure 1(b). Arrows in the 
figure show the directions of magnetic moments. The crystal structures of solid oxygen 
might be correlated with the magnetic moments even at high pressures. 

A Raman experiment reported a phase transition from a-0 2 to another monoclinic 
or orthorhombic structures below 3 GPa [2]. X-ray diffraction experiments using high- 
brilliance synchrotron radiation source, however, reported different phase diagrams. 
Akahama et al. reported a-02 transforms to e-02 directly at 7.2 GPa at 19 K [4]. 
Although they also found an anomaly in lattice constants just before the transition 
to e-C>2, it was not considered as a sign of a transition to other phases. Gorelli et al. 
also reported the stability of C1-O2 up to about 5.3 GPa, and they found orthorhombic 
5-0 2 (space group: Fmmm) above 5.3 GPa at 65 K [5|. Although the monoclinic a 
structure is ascertained to be stable up to about 5 GPa at low temperatures in both 
experiments, the stable structure at pressures between 5-7 GPa is considered to be still 
open question. For higher pressures, recent X-ray diffraction experiments have revealed 
that e-02 consists of Og clusters [HI [7j. The e phase transforms to the metallic ( phase 
at higher pressure (~ 100 GPa) and finally to a superconducting state [8j [9j [10]. The 
mechanism of the transition from or to the e phase and the structure of the ( phase are 
still unknown. 

Regarding theoretical studies, Etters et al. obtained the structure of a-02 using 
semi-empirical pair-potentials including magnetic interactions between O2 pairs and 
predicted a phase transition from the a to an orthorhombic phase at 2.3 GPa [3]. 
We performed theoretical calculations using ab-initio pair-potentials and reported 
that the monoclinic a structure is stable up to 6 GPa [11]. A part of the results 
will be presented in this paper. First-principles investigations based on density 
functional theory (DFT) [12] were performed especially for higher pressure phases e 
and ( [TU [HI [15j [TBI H3 HE]- The predicted structures for insulating e-02 is however 
not consistent with the experimental one. This is caused by a failure of the local 
density approximation (LDA) in describing the magnetic interaction between oxygen 
molecules JT9] . Because the failure does not influence the metallic state calculations, 
the structural relaxation based on DFT with an appropriate starting structure can 
predict the structure of the ( phase. A recent first-principles study reported a plausible 
structure for the metallic ( phase [IB] . In this paper, we study the insulating a and S 
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phases. The total energy is evaluated as a superposition of pair-potentials obtained by 
quantum chemistry calculations including configuration interactions. 

2. Theory and Calculation 

Under low pressures, solid oxygen consists of weakly bound oxygen molecules. Therefore 
the interactions between molecules need to be reproduced exactly in order to calculate 
the structure. The total energy calculation based on DFT is often used to investigate 
structures of various crystals. Unfortunately, DFT calculations within the LDA (or 
with gradient corrections) does not give correct intermolecular interactions for magnetic 
oxygen molecules. It is related with the symmetries of the highest occupied molecular 
orbital (HOMO) and lowest unoccupied molecular orbital (LUMO), and the well-known 
bandgap- underestimation problem of the LDA |19j . Both of the HOMO and LUMO 
of interacting oxygen molecules are originated in the 7r g (gerade) orbital of O2. Thus 
the underestimation of the HOMO-LUMO gap cause a serious overestimation of the 
exchange energy. We actually confirm DFT calculations do not reproduce the a 
structure at zero pressure. In many cases, quantum chemistry calculations including 
the configuration interactions reproduce intermolecular interactions correctly. In this 
paper, we evaluate the total energy as a superposition of intermolecular potentials (pair- 
potentials) calculated by an ab-initio method with configuration interactions. 

In order to reduce degrees of freedom of the O4 system, we introduce following 
assumptions into the calculations. 

(i) An antiferromagnetic ordering in the ab (basal) plane. 

(ii) The C2/m lattice symmetry. 

(iii) A molecular axis perpendicular to the ab plane. 

These assumptions are consistent with known properties of the structure of a-C>2, and 
we can treat the orthorhombic 5 structure as a special case that the monoclinic angle 
P* = 90°. 

In the orthorhombic structure, two types of the magnetic ordering along the c-axis 
are allowed: the ferromagnetic or anti-ferromagnetic orderings. In the a phase, the 
magnetic ordering along the c-axis is obviously ferromagnetic because of the periodicity 
along c-direction. Although recent neutron diffraction experiment reported anti- 
ferromagnetic ordering along the c-axis in the 5 phase [20] , we assume the ferromagnetic 
ordering along the c-axis because it simplify the treatment of the transition from a to 
5 phases. We have however confirmed that the difference in the total energy between 
the two kinds of orderings on the same lattice is very small. It is ascribed to the weak 
magnetic interaction between molecules in different c-layers. 

Following the assumption 1, there are two kinds of pairs of oxygen molecules. One is 
the pair in which magnetic moments are parallel (ferromagnetic pair, F). The other is the 
pair in which magnetic moments are anti-parallel (antiferromagnetic, AF). Assumptions 
2 and 3 restrict geometric configuration of molecules. In addition to these assumptions, 
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we fix the inter-atomic distance in O2. By these constraints, degrees of freedom of the 
O4 system, in other words the number of geometric parameters of pair-potentials, are 
reduced from twelve to two. Remaining degrees of freedom are the relative position of 
molecules. The relative position of oxygen molecules (x, z) is defined as given in figure 2. 
Consequently, the total energy is written as 

[/total = E P\xJ mn) Z n ) + U AF (xf£ n , Z n )\ . (1) 
l,m,n 

In the equation, pair-potentials U F and U AF denote interactions of the F and AF pairs, 
respectively. l,m and n are indices of lattice vectors and coordinates xf mn ,x A ^ n and z n 
are given as 

\J {la — nc cos f3) 2 + (mb) 2 , (2) 



x 



Iran 




2 / / ^\ \ 2 



4L = y(^ + |Ja-nccos/?J + {{ m +^) b ) » ( 3 ) 

z n = |ncsin (3\ . (4) 

Here, a, b, c and (3 are lattice constants and the monoclinic angle, respectively. Pair- 
potentials U F and U AF are calculated by a complete active space self-consistent 
field (CASSCF) method [21]. The CASSCF calculation is performed with GAMESS 
program [22] using 3s2pld atomic natural orbital basis set [23]. Details of calculations 
and obtained pair-potentials have been given in our previous brief report 



3. Results and Discussions 

3.1. Structures up to 6 GPa and its origin 

A Raman experiment reported a-0"2 transforms to another monoclinic structure or 
orthorhombic structure under 3 GPa [2]. X-ray diffraction experiments by Akahama et 
al. and Gorelli et al, however, showed the a structure is stable up to higher pressure. 
In this section, we present a result of the theoretical calculations that the a structure 
is stable up to 6 GPa. 

The calculated pressure dependences of lattice parameters up to 6 GPa are shown in 
figures 3(a), (b) and (c). The results of higher pressure, which is also shown in the figures, 
will be discussed later. In figures, solid lines present calculated lattice parameters, and 
closed circles, triangles and diamonds are experimental results [U H]. The values of 
pressure are evaluated by the numerical differentiation of the total energy with respect 
to the volume (P = —dE/dV). 

As previously mentioned, Raman and theoretical studies reported phase transitions 
in this pressure range [2j [3]. The structural transition predicted by the theoretical 
calculation was accompanied by abrupt changes of lattice parameters at 2.3 GPa. As 
shown in figure 3, present results decrease continuously with increasing pressure and 
there are no sign of the structural transition to the 5 structure up to 6 GPa. 

As shown in figure 3(a), the crystal is less compressible along the b- axis than the 
a-axis. The relation of the anisotropic lattice compressibility and the AF ordering 
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in the ah plane is discussed by Akahama et al [4j. The b/a ratio is 0.638 at zero 
pressure, and it increases with pressure and reaches 0.691 at 6 GPa. The experimental 
value is 0.635 at zero pressure p], and it reaches 0.687 at 6 GPa [I]. The agreement 
between the calculation and experiments is very satisfactory. As previously reported, 
the crystal structure of a-0 2 is considered as a result of magnetic interactions between 
oxygen molecules [3] . In order to clarify how the magnetic interactions affect the crystal 
structure, we consider a simple model of the ab plane as shown in figure 4. This figure 
shows the top view of the ab plane. The molecular axes are perpendicular to the figure, 
and arrows represent magnetic moment directions of O2 molecules. The unit cell contains 
ten pair of oxygen molecules: four AF pairs at a distance r = \/a 2 + b 2 /2 and three 
kinds of F pairs at distances a, b and 2r. For simplicity, we neglect interactions between 
molecules in which intermolecular distances are 2r. Consequently the (total) energy of 
this system is written as 



E(r, b) = AA(r) + 2 \F(^Ar 2 - b 2 ) + F(b)\ , (5) 

where A(x) and F(x) denote interactions (pair-potentials) of the AF and F pairs at a 
distance x. Calculated intermolecular potentials A(x) and F(x), which correspond to 
U AF (x, 0) and U F (x, 0) in (CQ), are shown in figure 5. Optimal energy of ([5]) are obtained 
at (a, b, r) = (5.42, 3.41, 3.20), and those are indicated by arrows in the figure. The 
obtained values are very close to the experimental values at zero pressure (5.40, 3.43, 
3.20) p. 

Evidently ([5]) is dominated by the first term because of its coefficient and the 
fact that A(x) has deeper minimum than F(x). Thus we obtain the optimal value 
of r, r = 3.2 where A(r) takes the lowest energy. Then the problem is simplified as 



E(r ,b) = AA(r ) + 2E 2 (b), where E 2 (b) = F(a(b))+F(b) and a(b) = ^JAr 2 ~b 2 . E(r ,b) 
takes a minimum on the condition dE ^ = 0, i.e., 

do 

a(6) dF|6) =6 dFW. ( 6 ) 
do da 

This condition is symmetric for b and a, then a = 6(= V2r ) satisfies the condition. It 
gives, however, the highest energy because is negative at b = \f 7 Lr§ ~ 4.5 as shown 
in figure 5. Another solution 6 of © is located at slightly larger side of the minimum 
of F(x): around 3.4. As a result, the optimal value of r and 6, the structure of the 
ab plane therefore, strongly depend on the position of the minimum of A(x) and F(x). 
The difference of A(x) and F(x) corresponds to the exchange energy. If there are no 
magnetic interactions between the molecules, A(x) is equivalent to F(x). In this case, 
following the above discussion, the optimal value of r is equal to 60 an d it gives the 
close-packed triangular lattice. The distortion is not introduced if there are only one 
kind of pair-potential. Namely the origin of the distorted triangular lattice in the a 
phase is attributed to the existence of the two kinds of pair-potentials, in other words, 
the magnetic interaction. 

The anisotropic lattice compressibility of the ab plane, which is pointed out by 
Akahama et al. jl], can be easily understood from figure 5. First we consider a situation 
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that the lattice is compressed along the a-axis and the length along the 6-axis is fixed. 
In this case, r decrease following the relation r = \/ a 2 + b 2 and the first term of (jSJ) 
increases by the repulsive interaction. Some part of the increased energy is however 
cancelled out by the term of F(a) because it is in the attractive region. On the other 
hand, when the lattice is compressed along the 6-axis and a-direction is fixed, both term 
of A(r) and F(b) increase, and as a result the total energy E increase steeply compared 
with the previous case. This is the origin of the anisotropic lattice compressibility of 
the a phase. 

3.2. Structures above 6 GPa 

X-ray diffraction experiments by Akahama et al. and Gorelli et al. showed the 
monoclinic a structure is stable up to higher pressure than the previously proposed 
transition pressure. Around 6 GPa, however, they proposed different structures. 
Akahama et al. reported a direct structural transformation from a- to e-0"2 at 7.2 
GPa. On the contrary, Gorelli et al. found 5-02 between a- and e-02- In this section, 
we discuss the stability of a- and 5-02- 

In figure 3, dashed lines above 6 GPa show lattice parameters of optimal structure. 
The dotted lines present extrapolated values from the data below 6 GPa, namely 
those form monoclinic a structure. Although optimal values (dashed lines) show the 
orthorhombic structure is stable above 6 GPa, the energy of the orthorhombic structure 
is very close to that of the monoclinic (dotted lines) structure. 

Since the <5-0 2 can be considered as a special case of a-0 2 (/3*=90°), we illustrate the 
total energy as a function of (3* and pressure. Figure 6 presents a pressure dependence 
of the cross section of the total energy from 4 GPa to 7 GPa. In the figure, lattice 
constants a, b and the volume of the unit cell are fixed to optimal values presented 
by the solid and dashed lines in figure 3(a). Therefore the independent parameter is 
only (3* at each pressure. At the low pressure range, the total energy takes minima at 
two values of (3*, which give equivalent monoclinic a structures, and the orthorhombic 
structure {(3* = 90°) is unstable. With adopted constraints of the crystal structure 
and geometric configurations of molecules, it can be derived analytically that a first 
derivative of the total energy with respect to (3* is zero at (3* = 90° (orthorhombic 
structure) at each pressure, i.e., the total energy takes a maximum or minimum values 
at the 5 structure [19J. With increasing pressure, the 6 structure becomes stable at 
around 6 GPa in place of the a structure, the energy difference between these structure 
is however less than 0.5 meV (~ 5 K) even at 7 GPa. Unfortunately, in the precision of 
present calculations, the energy difference may be too small to decide which structure is 
stable. We obtained however an important result that the monoclinic-angle-dependence 
of the total energy is very small in this pressure range. Akahama et al. and Gorelli 
et al. proposed different structures in this pressure range. It is reasonable because 
the experiments were performed at temperatures higher than the energy difference we 
obtained. The very small energy difference in this pressure range implies that several 
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structures which have various monoclinic angles can be allowed at finite temperature. 
Mita et al. performed Raman scattering experiments and observed peaks above 5 
GPa [21]. Although the origin of the peaks is unresolved, they indicated the existence of 
the complicated mixed phase at this pressure range. Further investigations are desired 
to determine the origin of the unknown peaks. 

At high temperatures, many experiments demonstrated the stability of 
orthorhombic 5-02- The stability of 5-02 at high temperature is understood from 
figure 6. As shown in the figure, the energy surface forms a double-well or parabolic 
shape. In either case, they are symmetric with respect to (3* around (3* = 90° and the 
energy barrier separating the wells are very low. Therefore the expectation value of f3* 
becomes 90° even at low pressure because of the effect of thermal excitations. 

Akahama et al. reported abrupt changes of lattice parameters at the phase 
boundary of the a and e phases. In the present calculations, however, corresponding 
change of lattice parameters are not obtained up to 12 GPa. This may be caused by 
constraints on the crystal and/or magnetic structure in our calculations. Recent neutron 
diffraction experiment shows a disappearance of the magnetic long-range ordering in the 
e phase, but it does not exclude the short-range order [25]. Intermolecular potentials 
adaptable to other magnetic configurations may be required to obtain the e structure. 

4. Summary 

We investigated crystal structures of solid oxygen under pressure with ab-initio pair- 
potentials taking account of the difference of spin states. Under 6 GPa, the obtained 
structure is monoclinic a and calculated lattice parameters agree well with the results 
of the X-ray diffraction experiment. The origin of the distortion of the basal plane and 
that of the anisotropic lattice compressibility of a-02 are clearly demonstrated that they 
are caused by the existence of two kinds of intermolecular interaction, namely magnetic 
interactions. The crystal structure under the pressure range from 6 to 7 GPa is discussed 
using the total energy surface projected on the P-/3* space. Although orthorhombic 5- 
O2 is obtained as the stable structure, the energy difference between a- and 5-O2 is less 
than 5 K in this pressure range. This implies a possibility of a mixed phase including 
some structures which have various monoclinic angles. It seems to be consistent with 
recent X-ray diffraction and Raman experiments. 
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FIGURE CAPTIONS 



figure 1 (a) Monoclinic unit cell of CK-O2. The unit cell includes two oxygen 
molecules. At (3* = 90°, the unit cell becomes an orthorhombic Fmmm structure, 
(b) The basal plane of a-02- Oxygen molecules are denoted by circles and molecular 
axes are perpendicular to the ab plane. Arrows represent magnetic moment 
directions. 

figure 2 A schematic view of coordinates for pair-potentials. The vector (x, z) is 
the relative position of molecules. 

figure 3 Pressure dependences of (a) lattice constants and monoclinic angles (b) 
(3 and (c) (3* up to 7 GPa. Calculated values are denoted by solid lines and 
experimental values are represented by closed circles, triangles and diamonds. See 
text regarding the dashed and dotted lines. 

figure 4 A model of the ab plane. Circles denote oxygen molecules, and arrows 
represent magnetic moment directions. The unit cell includes four AF pairs and 
six F pairs. Two F pairs in which intermolecular distances are 2r is not taken into 
account in the discussion. 

figure 5 The optimized values of (a, b, r) and pair-potentials. These pair-potentials 
correspond to U F (x, 0) and U AF (x, 0), which are defined in previous section. Arrows 
denote the optimal values of (a, b, r). Note that the values of b and r are at almost 
minimum of pair-potentials. 

figure 6 The contour plot of the total energy as a function of the pressure and (3* 
in the pressure range 4-7 GPa. The volume of the unit cell and lattice constants a 
and b are fixed at optimal values at each pressure, which are denoted by solid and 
dashed lines in figure 3. The difference in the energy between the monoclinic a and 
orthorhombic 5 structures is very small at around 6 GPa. Contours increase by 2 
meV. 




g 

..(J! 

9 6 

(0, 0) ^ 

O Z V 



(x, z) 



Figure 2. 




Figure 3. 





2 2.5 3 3.5 4 4.5 5 5.5 6 
INTERMOLECULAR DISTANCE [A] 



Figure 5. 



110 




PRESSURE [GPa] 



Figure 6. 



